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Abstract 

The scope of the present study is Eulerian modehng and simulation of polydisperse liquid sprays 
undergoing droplet coalescence and evaporation. The fundamental mathematical description is the 
Williams spray equation governing the joint number density function /(w,u;x, t) of droplet volume 
and velocity. Eulerian multi-fluid models have already been rigorously derived from this equation in 
Laurent et al. [22]. The first key feature of the paper is the application of direct quadrature method of 
moments (DQMOM) introduced by Marchisio and Fox [24] to the Williams spray equation. Both the 
multi-fluid method and DQMOM yield systems of Eulerian conservation equations with complicated 
interaction terms representing coalescence. In order to focus on the difficulties associated with treat- 
ing size-dependent coalescence and to avoid numerical uncertainty issues associated with two-way 
coupling, only one-way coupling between the droplets and a given gas velocity field is considered. 
In order to validate and compare these approaches, the chosen configuration is a self-similar 2D ax- 
isymmetrical decelerating nozzle with sprays having various size distributions, ranging from smooth 
ones up to Dirac delta functions. The second key feature of the paper is a thorough comparison of the 
two approaches for various test-cases to a reference solution obtained through a classical stochastic 
Lagrangian solver. Both Eulerian models prove to describe adequately spray coalescence and yield 
a very interesting alternative to the Lagrangian solver. The third key point of the study is a detailed 
description of the limitations associated with each method, thus giving criteria for their use as well as 
for their respective efficiency. 
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1 Introduction 



In many industrial combustion applications such as Diesel engines, fuel is stocked in con- 
densed form and burned as a dispersed liquid phase carried by a gaseous flow. Two phase 
effects as well as the polydisperse character of the droplet size distribution (since the droplets 
dynamics depend on their inertia and are conditioned by size) can significantly influence 
flame structure. Size distribution effects are also encountered in a crucial way in solid pro- 
pellant rocket boosters, where the cloud of alumina particles experiences coalescence and 
become polydisperse in size, thus determining their global dynamical behavior [17,18]. The 
coupling of dynamics, conditioned on particle size, with coalescence or aggregation as well 
as with evaporation can also be found in the study of fluidized beds [36] and planet formation 
in solar nebulae [5,6]. Consequently, it is important to have reliable models and numerical 
methods in order to be able to describe precisely the physics of two-phase flows where the 
dispersed phase is constituted of a cloud of particles of various sizes that can evaporate, coa- 
lesce or aggregate and also have their own inertia and size-conditioned dynamics. Since our 
main area of interest is combustion, we will work with sprays throughout this paper, keeping 
in mind the broad application fields related to this study. 

Generally speaking, two approaches for treating liquid sprays corresponding to two levels 
of description can be distinguished. The first, associated with a full direct numerical simu- 
lation of the process, provides a model for the dynamics of the interface between the gas 
and liquid, as well as the exchanges of heat and mass between the two phases using vari- 
ous techniques such as the Volume Of Fluids (VOF) or Level Set methods [3,15,19,35]. This 
"microscopic" point of view is very rich in information on the detailed properties at a more 
local level concerning, for example, the resulting drag exerted on one droplet depending on 
its surroundings. The second approach, based on a more global point of view, describes the 
droplets as a cloud of point particles for which the exchanges of mass, momentum and heat 
are described globally, using eventually correlations, and the details of the interface behavior, 
angular momentum of droplets, detailed internal temperature distribution inside the droplet, 
etc., are not predicted. Instead, a finite set of global properties such as mass, momentum, 
temperature are modeled. Because it is the only one for which numerical simulations at the 
scale of a combustion chamber or in a free jet can be conducted, this "mesoscopic" point of 
view will be adopted in the present paper. 

Furthermore, we are interested in sprays where droplet interactions (e.g., coalescence) have to 
be taken into account, which corresponds to liquid volume fractions between 0.1% and 1%. 
O'Rourke [30] classified the various regimes from the "very thin spray", which are trans- 
ported by the gaseous carrier phase without influencing the gaseous phase, through the "thin 
spray" regime, for which there is two-way coupling through the momentum equation be- 
tween the two phases, up to the "thick spray" regime for which the volume fraction of liquid 
is high enough so that droplet-droplet interactions have to be taken into account, but is still 
low enough so that the liquid volume fraction is negligible as compared to the gaseous one. 
Because our primary focus is on the ability of Eulerian methods to capture droplet coales- 
cence, our study is limited here to the "thick spray" regime. By restricting our attention to 
one-way coupling, we can avoid difficulties (e.g., grid convergence) associated with using 
Lagrangian methods with two-way coupling, and it will thus be possible to make detailed 
comparisons between Eulerian and Lagrangian simulation results. 
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In the mesoscopic framework, there exists considerable interest in the development of numer- 
ical methods for simulating sprays [18,17,27,28,22,32]. The principal physical processes that 
must be accounted for are (1) transport in real space, (2) droplet evaporation, (3) acceleration 
of droplets due to drag, and (4) coalescence of droplets leading to polydispersity. The major 
challenge in numerical simulations is to account for the strong coupling between these pro- 
cesses. Williams [37] proposed a relatively simple transport equation based on kinetic theory 
that has proven to be a useful starting point for testing novel numerical methods for treating 
coalescing liquid sprays. In the context of one-way coupling, the Lagrangian Monte-Carlo 
approach [9], called Direct Simulation Monte-Carlo method (DSMC) by Bird [4], is gener- 
ally considered to be more accurate than Eulerian methods for solving Williams equation. 
However, its computational cost is high, especially in unsteady configurations. Moreover, in 
applications with two-way coupling, Lagrangian methods are difficult to couple accurately 
with Eulerian descriptions of the gas phase. There is thus considerable impetus to develop 
Eulerian methods for describing sprays. In this paper, we limit our attention to one-way 
coupling with a given (laminar) gas velocity field (i.e., one-way coupling with a given gas 
velocity field.) Thus no turbulence models are required to close the spray equation. 

In a recent paper Laurent et al. [22] have demonstrated the capability of an Eulerian multi- 
fluid model to capture the physics of polydisperse evaporating sprays with one-way coupling. 
This approach relies on the derivation of a semi-kinetic model from the Williams equation us- 
ing a moment method for velocity, but keeping the continuous size distribution function. This 
distribution function is then discretized using a "finite-volume" approach that yields conser- 
vation equations for mass, momentum (and eventually other properties such as temperature) 
of droplets in fixed size intervals called "sections" extending the original work of Tambour, 
Greenberg and collaborators [12,13]. Even though this approach has recently been extended 
to higher order by Laurent [20] and Dufour [7,8], the necessity to discretize the size phase 
space can be a stumbling block in some applications. Moment methods, on the other hand, 
do not encounter this limitation. 

In this work, we apply the recently developed direct quadrature method of moments (DQ- 
MOM) [24] to treat Williams equation in a Eulerian framework. As its name implies, DQ- 
MOM is a moment method that closes the non-linear terms (e.g., droplet coalescence) using 
weighted quadrature points (abscissas) in phase space. Such a closure relates to the construc- 
tion of an approximated number density function from a set of moments under the form of a 
sum of Dirac delta functions, the support of which corresponds to the abscissas. However, it is 
important to make a clear difference between such an Eulerian approach and the correspond- 
ing Lagrangian approach, for which the number density is approximated by a large number of 
numerical "parcels". The evolution of abscissas and the corresponding weights are governed 
by the dynamics of a few moments, whereas the evolution of the parcels are governed by the 
Williams equation since they are a stochastic discretization of this equation. Consequently, 
the DQMOM usually involves a very restricted number of unknowns on a Eulerian mesh, 
whereas the Lagrangian method involves a very large number of unknowns that are followed 
along their trajectories in phase space. 

The DQMOM method distinguishes itself from other quadrature methods (e.g., QMOM 
[26,25]) by solving transport equations for the weights and abscissas directly (instead of 
transport equations for the moments). The source terms for the transport equations depend on 
the physical processes involved. For Williams equation, we show in Section 2 that laminar 
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transport and drag result in source terms that are independent of the choice of moments and, 
in fact, are equivalent to those used in Lagrangian formulations. When evaporation does not 
lead to the disappearance of droplets in finite time, this is also true for the evaporation pro- 
cess. On the other hand, coalescence leads to a linear system for the source terms for which 
the coefficient matrix depends on the choice of moments. The applicability of DQMOM to 
Williams equation thus depends on whether or not a particular choice(s) of moments can be 
found that leads to a non-singular linear system. When the evaporation law allows the dis- 
appearance of droplets in finite time the equations for the moments of the number density 
function not only involve unclosed integral terms, but also the flux of disappearing droplets, 
i.e. the pointwise value of the number density function at zero size. This quantity has then to 
be closed since it has a strong influence on the dynamics of the whole set of moments; it leads 
to a significant difficulty since it corresponds to the reconstruction of a pointwise value of the 
number density function from a set of its moments. In this study, we propose a solution to 
this difficult issue. Note that because spatial transport is treated explicitly, it suffices to tackle 
the flux problem in the homogeneous case. We will see that a key point is to provide a flux 
closure that yields stable moment dynamics and a non-singular linear system in the DQMOM 
framework. 

Let us also underline that the transport terms in the systems of conservation equations for 
both Eulerian models are the same and given by pressureless gas dynamics. The structure 
of these transport terms and the associated difficulties have been the subject of several stud- 
ies and there are numerical methods designed in order to treat the resulting singularities as 
shown in [22]. The question of the computational efficiency of such Eulerian approaches (es- 
pecially in coalescing systems) is a key question since these methods are intended to be used 
in more realistic unsteady configurations as an alternative to the too costly Lagrangian meth- 
ods for poly disperse sprays. We have already studied this question in [22] where the Eulerian 
multi-fluid approach was shown to offer a good precision with a relative low cost [22]. Be- 
cause of the similarity of the transport terms for both Eulerian approaches, the conclusions 
about the computational efficiency presented in [22] are also valid for the DQMOM method. 
Consequently we focus our study and comparisons on stationary configurations for which 
we are sure to have a reference solution at our disposal and from which we can obtain firm 
conclusions about the capabilities of the various approaches. 

In Section 3, we present the chosen test configuration, which is a self-similar 2D axisym- 
metrical decelerating nozzle and sprays with two inlet distributions: a smooth monomodal 
function and Dirac delta functions. We also discuss in detail the reasons (e.g., significant co- 
alescence rates) for the choice of the test cases, and why they are particularly challenging 
for the various numerical methods. Finally the Lagrangian solver, the numerical subtleties 
for obtaining the associated reference solution, as well as the multi-fluid method are then 
presented. In Section 4, we consider the results for the various test cases including combina- 
tions of coalescence, linear evaporation in terms of volume (since it conserves the number of 
droplets and thereby eliminates the need to model the evaporative flux) and the usual non- 
linear evaporation law (for which the evaporative flux must be modeled.) We present results 
for the most difficult test cases, designed to highlight the challenges one would encounter in 
more realistic cases. The results are compared to a reference solution obtained through a La- 
grangian stochastic algorithm [17]. The advantages and limitations of the Eulerian methods 
are then analyzed in detail in terms of precision and efficiency. It is shown that the DQMOM 
method offers very interesting features in a number of situations (e.g., strongly coalescing 
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droplets), and is a good candidate for more complex configurations. 



2 DQMOM for Williams equation 



The Williams transport equation [37] for the joint volume, velocity number density function 

f{v,u; X, t) is 

dtf + u- d^f + a (Rvf) + 9u • (F/) = r, (1) 

where is the evaporation rate, F is the drag force acting on the droplet, and T is the 
coalescence term. Note that specific forms for the evaporation rate and drag law are not 
required for DQMOM. However, in this work we will consider one-way coupling with a 
given gas velocity that appears in F. Using standard assumptions [22], we can write the 
coalescence term in two parts: F = + Q'^u where 



Q:oii = - 11 Bi\u-u*\,v,v*)fiv,u)fiv*,u*)dv*du\ (2) 
Q^u [ B{\n'> - n*\y, v*)f{v\ n^)f{v*, n*)Jdv* du*, (3) 

= V — V*, u* = {vu — v*u*)/ {v — V*), and J = (v/v^)^ is the Jacobian of the transform 
{v, u) (f u*) with fixed (f *, u*). The collision frequency function B is defined by 

5(|u-u*|,^;,^;*) = £;eoai(|u-u*|,^;,^;*)/3(^;,^;*)|u-u*|, (4) 

where E'coai is the coalescence efficiency probability, which, based upon the size of droplets 
and the relative velocity, discriminates between rebound and coalescence, and 



(3{v, V*) — TT 



\47r/ V 47r y 



(5) 



For simplicity, we will take £^coai = (no coalescence) or Ecoai = 1; however, any other 
functional form could be used in the derivation that follows. A more general version of the 
spray equation would include the droplet temperature and molecular composition. For sim- 
plicity, we consider only the volume and velocity in this work. Finally note that adding spatial 
diffusion terms in Eq. (1) would generate additional terms in DQMOM [24]. 

One of the principal mathematical difficulties when developing Eulerian solvers for Eq. (1) 
is the accurate treatment of the coalescence term. Indeed, the integral form of T leads to 
highly non-local and non-linear interactions in volume-velocity phase space. A "direct" Eu- 
lerian solver would require discretization of the high-dimensional phase space (in addition 
to real space), and would thus be computationally intractable. In contrast, multi-fluid models 
discretize only the volume phase space and use the average velocity conditioned on droplet 
size (i.e., the mono-kinetic assumption [21]), while moment methods (such as DQMOM) 
provide closures based on a finite set of moments. Before applying DQMOM to Eq. (1), we 
should note that the coalescence term is defined such that the moments representing mass and 
momentum are conserved: 

J pvr{v,u)dvdu^O (6) 
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and 



j pvur{v, u) dv du = 0, (7) 

where the liquid density p is assumed to be constant. These conservation properties must be 
retained in numerical approximations used to treat Eq. (1) (as we shall see is the case with 
DQMOM). 

The DQMOM approximates the density function by weighted delta functions in volume- 
velocity phase space [11,24]: 

N 

f(v,u) ^^WnS{v -Vn)S(u-Un) (8) 

n=l 

where 6{u — u„) = 5{ui — Mi,„)5(m2 — M2,n)^('"3 — ""s.n)- Note that in this formulation, 
the weights w„ and abscissas (Vn, u„) are Eulerian fields. Application of DQMOM results 
in closed transport equations for the number density, mass density, and momentum density, 
respectively, of the form: 

dtWn + 9x ■ (w„u„) = a„, (9) 

dt (WnPVn) + 5x • {WnPVnnn) = pbn, (10) 

and 

dt (WnPVnnn) + 8^ ■ (w„/W„U„U„) = pC„, (11) 

where a„, and Cn are source terms that are found from the right-hand side of Eq. (1) 
as described below. These equations can be solved with appropriate initial and boundary 
conditions to find the fields w„(x, t) and (f„(x, t), u„(x, t)) appearing in Eq. (8). Note that 
Eqs. (9-11) are equivalent to an Eulerian multi-fluid model [22], but with the source terms 
on the right-hand side determined using DQMOM. 

The DQMOM approximation for the moments of the number density function are found 
directly from Eq. (8): 

N 

{v^u{u^ul) = / v%u^ulf{v, u) d^; du = ^ <„. (12) 

n=l 

The fundamental idea behind DQMOM is that we should choose the weights and abscissas 
such that as many moments as possible are determined by the moment transport equations 
found from Eq. (1). Note that there are a total of N weights, volume abscissas, and 3A^ 
velocity abscissas and (equivalently) 5A^ unknown source terms in Eqs. (9-11). We will thus 
need to choose 5A^ independent moments to determine the source terms. We will return to the 
subject of how to choose the moments in Section 2.4. The procedure for using these moments 
to find the source terms is described next. 
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2. 1 Space and time derivatives 



The space and time derivatives in Eq. (1) generate the corresponding terms in Eqs. (9-11). 
These are found by formally inserting Eq. (8), and differentiating: 

N 

dtf + U- (9x/ ='^S{V - Vn)S{u - Un) [dtWn + 0^ ■ (UnW„)] 
n=l 

N 

- WnS''^\v - Vn)5{u - U„) [dtVn + U„ ■ d^Vn] 
n=l 

N 

- WnS{v - Vn)S^^\u - U„) ■ [dtUn + U„ • d^Un] (13) 
n=l 

where (5*^^^ (ip) = d6{ilj) / dip and 5^^^ {ip>) is a vector with components Si'^ (ip) = 5*^^^ (■?/'i)5(-?/'2)^ ('V's); 
Si^\il)) = 6{ilJi)6^^\ilj2)6{'ilJs), and ^^HV') = 5(t/'i)5(^/'2)5(^H^/'3). Using the definitions of 
the source terms, Eq. (13) can be rewritten as 



dtf + u- d^f = Y,[S{v - Vn)6{u - u„) + VnS'^^\v - Vn)6{u - u„) a„ 

n=l 

N 

- J2 [S^^\v - Vn)S{n - Un) - V:^^6{V - Vn)S^^\u - U„) • Un 6n 
n=l 

N 

-J2'''n'S{v-Vr,)S^'\n-U^)-Cr,. (14) 



n=l 

Note that this expression is linear in the source terms (a^, 6„, c^). 

The next step is to apply the moment transform to Eq. (14). Formally, this yields 

N 

/ v%u^ul {dtf + u • d^f) dvdu = i:(l - A;)^;X,n<n<n«« 

n=l 

N 

N 



n=l 

N 



n=l 



where, unless otherwise noted, the definite integrals cover all of phase space. The next step is 
to consider the terms in Eq. (1) that correspond to transport in volume-velocity phase space. 



2.2 Phase-space transport 

We begin by rewriting Eq. (1) as 

dtf + u-d^f^P, (16) 
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where the phase-space transport terms are defined by 



P=-d,{R,f)-d^-{Ff) + r. (17) 
We can then define the moment transform of the phase-space terms by 

P{k,l,m,p) = J v''u[u'^ulP dvdu. (18) 

Note that if the moments P{k, I, m, p) are known, Eq. (15) forms a linear system that can be 
solved to find the unknown source terms. We can compute the phase-space moments using 
Eq. (17): 

P{k, l,m,p) = -J v''u[u"^'ul [d, {RJ) + ■ (F/) - r] d^; du. (19) 

As shown next, the integrals on the right-hand side can be expressed in terms of the weights 
and abscissas, and a flux term corresponding to disappearance of droplets due to evaporation. 

Starting with the evaporation term in Eq. (19), we can use integration by parts to find 

/ v% (RJ) dv = -SkoRviO: u)/(0, u) - / kv''-^R,{v, u)fdv, (20) 
Jo Jo 

where Sko is the Kronecker delta. Using Eq. (8) in the final integral, we find 

I v^u'^uld, (RJ) dvdu = -Sko^{t)u'f^uf^u% 

N 

- kWnV^-'u[^^U^^ys^MVn, U„), (21) 
n=l 

where V'(^) is the evaporative flux of droplets at zero size and u/ is the velocity of droplets 
with zero volume (which will normally correspond to the fluid velocity). Note that the first 
term on the right-hand side of Eq. (21) will be non-zero only for k = 0, and corresponds to 
the loss of droplets due to evaporation. A fundamental question when applying DQMOM to 
evaporation problems is how to determine 'ijj{t) from the weights and abscissas. The value of 
'ip{t) corresponds to the value of the number density function at zero size, and in the case of 
the (P evaporation law, it is precisely the value of the number density as a function of droplet 
surface, which has no reason to be zero in general. Determining the value of ip{t), a pointwise 
information, from the values of moments is clearly a difficult task, for which we will propose 
a solution in the next subsection. On the other hand, the second term on the right-hand side 
of Eq. (21) is non-zero only for k > 0, and appears in closed form. 

Turning next to the drag-force term in Eq. (19), we can use integration by parts to find 

J u%. {Fjf) du^-J yr^F^f du for j = 1, 2, 3. (22) 
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Thus, the drag-force term becomes 



j v''u[u^uld^ ■ (F/) d^;du = 



N 



n=l 



^«l,n^l(^n, U„) + mU^]^F2{Vn, U„) + pU^l^F^.iVn, U^)] • (23) 



Note that this term appears in closed form. 

Turning now to the coalescence term, we will treat each of the two parts Q~oii and Q'loii 
separately. The first part yields in a straightforward manner 



/ 



N N 



v^v^^u'^ulQ^ii dvdu^-Y^Y. '^nWgV^ui^^u'^^y^^^B (I u„ - u, I , ^;„, ^;,) . (24) 



n=l q=l 



The second part requires a change in the order of integration, and a change of variables: 



J h{v,u)B{\u^ -u*\,v^,v*)f{v'',u'')f{v*,u*)Jdv*^ dvdu*du 
= J I^J^ h{v,u)B{\u'' -u*\,v'',v*)f{v^,u'')f{v*,u*)Jdv^ d^;*du*du 

J \ V* + V* / 

S(|u* - u*|, V*, v*)f{v^, u^)f{v*, u*) dv* dv^'du* du^ (25) 



where /i is an arbitrary function of v and u. It then follows that 



^ N N / 

n=l q=l \ 



Vn + Vg 



VnU2,n + VgU2,q \ ( VnUs^n + ^g^S.q 



Vn + Vq 



Vn + Vq 



B{\\ln -Ug\, Vn,Vq). (26) 



Note that the right-hand side of this expression is in closed form. 

Collecting together all of the terms, the moments appearing on the right-hand sides of Eqs. (21- 
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26) become 



TV 



n=l 



+ E 



n=l 



N N 



^ n=l 0=1 



Vn + Vq 



5(|U„ - Ug|,V„,Vq). (27) 



Note that due to the form of the coalescence term, the moments conserve mass (P(l, 0, 0, 0) = 
0) and momentum (P(l, 1, 0, 0) = P(l, 0, 1, 0) = P(l, 0, 0, 1) = 0) when evaporation and 
drag are null. Thus, the weights and abscissas in the DQMOM representation will keep the 
same conservation properties as the original model (i.e., as Eq. (1)). 

Comparing the terms in Eqs. (15) and (27), we can note that the evaporation and drag terms 
in the DQMOM representation can be solved for explicitly. Thus, the source terms can be 
written as 



bn = bl + WnRv{Vn,Un), 

Cn = < + WnUnRv{Vn, U„) + W„U„F(v„, U„), 



(28) 
(29) 



where source terms a^, 6* and c* in the transport equations are found by solving the linear 
system 



N 



N 



n=l n=l 

N 

+ j:<~'<n<,n<n (K>I,n + ^<4n + P^n) = P*ik,l,m,p), (30) 

n=l 

with the right-hand side given by 



N N 



P*{k, I, m,p) = -5k0lljU^flUj2U% + ^ E E "^nWq 



{Vn + Vq) 



k / VnUl^n + VqUi^q \ / VnU2,n + VqU2,q \ / VnU^^n + VgU^^q \ 



Vn + Vq J \ Vn + Vq 



VnV'l,nV%nV'3,n VgU[gU^gU^g 



Vn + Vq J 

B(\Un-Uq\,Vn,Vq). (31) 



The expression for the source terms (Eq. 30) completes the derivation of the DQMOM trans- 
port equations for the Williams spray equation. 



10 



In the absence of coalescence, Eq. (31) is particularly simple. Thus, the pure evaporation 
case is an interesting limit case for which a„, 6* , and c* will be non-zero only if the evap- 
orative flux 'ip is non-zero. However, the evaporative flux cannot be determined by moment 
constraints alone (see Section 2.3). If the evaporative flux is assumed to be null, the zero- 
order moment will remain unchanged in the absence of coalescence as long as some abscissa 
crosses the zero size limit and yields a pointwise singular and infinite flux as in Lagrangian 
methods when some parcels reach the zero size limit. However, as mentioned in the Introduc- 
tion, since there are only a few abscissas that describe the moment dynamics, such a singular 
behavior is not ideal for smooth number density functions (whereas it is the correct one if the 
number density function is a sum of Dirac delta function from the beginning as in the bimodal 
case that will be studied later). Consequently we need an evaluation of this flux function that 
guarantees a smooth flux as a function of time for smooth distribution functions. Even when 
coalescence is included, the moments may be poorly estimated if the evaporative flux is ne- 
glected. An example of such behavior can be found in the work of Mossa [29] where the 
droplet size distribution is presumed to be log-normal and where the evaporative flux at zero 
size is neglected, leading to numerical difficulties and a poor prediction of the second mo- 
ment. Thus, we will use a separate procedure, described next, to approximate the contribution 
due to the evaporative flux that yields a continuous in time flux, as well as a guarantee that 
the abscissas never cross the zero size limit. 



2.3 Evaporative flux 



The source terms cannot be computed directly from the moment constraints in Eq. (31) be- 
cause the evaporative flux is unknown. We must therefore apply additional (or different) con- 
straints to determine all of the unknowns. Considering only evaporation and setting drag and 
coalescence to zero in the right-hand side of Eq. (30), we obtain the following linear system: 

N N 

n=l n=l 

N 

n=l 

+ Skou'f,uf,u%^ ^ (32) 

with 5A^ + 1 unknowns a„, 6* , c* and ip. Note that because the right-hand side is null, only 
trivial solutions can be found using moment constraints. We will therefore introduce ratio 
constraints of the form 

^(^] =0, ^(^] =0 and ^(^] =0, 

Dt V^n+l/evap V^n+iyevap \Ujn+l J ^^^^ 

which are applied only for the changes due to evaporation. These constraints are motivated 
by the behavior of the weights and abscissas corresponding to sufficiently smooth and con- 
tinuous density functions. For example, if the surface density function is exponential and the 
evaporation rate is proportional to the surface area of a droplet, then the abscissas remain con- 
stant and the weights decrease monotonely. On the other hand, for singular density functions 
(e.g., composed of delta functions), the ratio constraints are expected to perform poorly. We 
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will look more closely at this issue in Section 4. 



It can be observed that the choice ofk = l = m= p = Om Eq. (32) leads to 

TV 

ilj = - ^ttn- (33) 

n=l 

Thus, the evaporative flux depends only on a„. Note that physically t/' > 0. Hence, if the value 
computed for t/j from Eq. (33) is negative (which is possible for very general evaporation 
rates), then a„, 6* , c* and ^jJ are set equal to zero. However, for the evaporation rate considered 
in this paper, it can be shown that the flux will be non-negative. 

Conservation of mass (A; = 1 and I — m — p — Oin Eq. (32)) leads to 

E^n = 0. (34) 

n 

Applying the ratio constraint for the abscissas yields 

Wn+iVn+ib*^ - WnVnbl+i = K forn = 1, . . . , - 1; (35) 
where the right-hand side is defined by 

En = WnWn+1 [VnRv{Vn+l) " Vn+lRv{v„)] ■ (36) 

Note that in order for there to be an evaporative flux, we will normally have En > for 
all n (assuming that vi < V2 < ■ ■ ■ < vn)- The case where En — occurs when Rv{v) is 
proportional to —v (i.e., the evaporation rate is proportional to the droplet volume). The more 
common case where En > occurs when Ry{v) is proportional to —v^^^ (i.e., the droplet 
surface area decreases linearly). In general, Rv{v) oc —v"' with 7 < 1 will lead to positive 
En- The physical interpretation for this difference is that for 7 < 1 the droplets will disappear 
due to evaporation in a finite time, while for 7 > 1 the disappearance time is infinite. The 
linear system formed from Eqs. (34) and (35) can be solved separately to find 6* . 

Conservation of momentum (k — 1 and m, or p = 1 in Eq. (32)) leads to 

TV 

E < = 0- (37) 

n=l 

Likewise, the ratio constraint for each component of the velocity yields 

Wn+lVn+lUjn+lC*n " WnVnUjnC*n+l = UjnUjn+lEn foV U = 1, . . . , N - 1. (38) 

Together with Eq. (37), this equation can be solved separately for each component (j — 
1, 2, 3) to find c;. 

The ratio constraint for the weights yields N — 1 equations for a^: 

Wn+lttn - Wnttn+l ^ fOT 71 ^ 1, . . . , N - 1. (39) 

Note that this constraint is satisfied by a„ = aWn where a is unknown. We must therefore 
choose one independent moment in Eq. (32) in order to solve for a. Since 6* and c* are 
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known, we can rearrange Eq. (32) as 



N 

n=l 
N 



Wr, 



T.(k - I - m - p)vt'u[,n<,ry3,nb: 
n=l 

N 

+ f^n~^'>J'in<,ntJ'in {}u^^nCl,n + ^t^2,li4,n + PU3^n4,n) > (40) 
n=l 

which can be solved with A; 7^ 1 to find a. If we choose, for example, k — 2 and I — m — 
p = as the independent moment, then the constraint becomes 

N IN 

a = 2 ^ / ^ ^;>„ (41) 

n=l / n=l 

and a depends only on 6*. However, if we choose k = 2 and I = m = p = 1, then the 
constraint becomes 

N 



n=l 

/ N 



vluinU2nU3nWn ■ (42) 



n=l 



For this choice, a is independent of u/. A choice that leads to a fully coupled system isk — 2, 
I — 2, m — p — 0, which yields 

N IN 

« = 2 ^ VnUinCln / '"l'^l,n'^n (43) 
n=l / n=l 



or/c = m= p = and I — 1, which yields 

N IN 



« = ^ (^*l,n&n " ^1,^) / II («l,n " Uf^) ■ (44) 



n=l / n=l 



Note that when w„ — > 0, we have ■Ui „ ^ Ufi and ^ ^fiK^ hence, this last constraint 
is consistent with this limiting behavior. These choices are asymmetric in the velocity com- 
ponents, and thus do not treat all components the same. A "symmetric" choice with similar 
properties is A; = 2 and I — m — p — 2ork — and I — m — p — 1, which lead to a 
more complicated constraint. The "best" choice will most likely be problem dependent. In 
our test cases, the choices with k = 2 give similar results, better than the ones with k = 0. 
The calculations are thus done with the value of a given in Eq. (41): this value is the simplest 
and can be shown to be non-positive as soon as > 0, at least for the case N — 2. 

In summary, the contribution due to evaporation is estimated by first solving separate linear 
systems for 6* and c* . The estimate for a„ = aWn is found using an independent moment 
constraint from Eq. (40) to find a. Finally, the evaporative flux ip is computed from Eq. (33), 
and should be non-negative. If ^ is negative (or equivalently if a is positive), then the con- 
tribution due to evaporation is null. The contribution due to coalescence is found by solving 
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a linear system of the form of Eq. (30) where the right-hand side is given by Eq. (31) with 
il) = 0. As described below, the final source terms (a„, 6* , c*) are found simply by adding 
together the contributions from the evaporative flux and coalescence. 

2. 4 DQMOM linear system 

The DQMOM representation of Williams spray equation is given by the transport equations 
for the weights and abscissas (Eqs. (9-1 1)). The source terms for these equations are found by 
solving the linear system as described above. The exact form of the linear system depends on 
the choice of moments. This choice, in turn, will determine if the system is well defined in the 
sense that the coefficient matrix is non-singular. A choice of moments that is consistent with 
the mono-kinetic assumption used in the multi-fluid model is to consider only moments of 
orders zero and one in the velocity components (i.e., l,m,p e {0, 1}). In this work, in order to 
make direct comparisons with the multi-fluid model, we will limit our consideration to such 
moments. In general, this choice of moments should allow for the best possible description 
of droplet coalescence, while at the same time ensuring that droplet mass and momentum are 
conserved. 

A choice of 5A^ moments that has been found to be non-singular is 





-l)/3 


ie{l,...,2N} 


with 


I — m — p — 






k — i 


ie{l, 


..,N} 


with 


I — 1, m — p 


= 




k — i 


ie{l, 


..,N} 


with 


m — 1, I — p 


= 




k — i 


ie{l, 


..,N} 


with 


p — 1, I — m 


= 0. 


(45) 



For N > 2, this choice of moments includes the surface area and the volume of the droplets, 
which are important variables for evaporating spray, as well as their momentum. The linear 
system can then be written in matrix form (showing only non-zero components) as 



Ai 


A2 


El 


E2 


E3 




a 




Pa 


A3 


A4 










b* 




Pb 


Bi 


Ci 


Di 












Pi 


B2 


C2 




D2 






c* 




P2 


B3 


C3 






D3 




C3 




P3 



where the matrices Aj, Bj, Cj, Dj and are all NxN, and a, b* and c* are column vectors 
formed from the components a„, 6* and c*„, respectively. In general, the exact definitions of 
the other matrices will depend on which constraints are used to define the system, i.e., Eq. (30) 
or those described in Section 2.3. Nevertheless, the form of the linear system is the same in 
all cases. As noted earlier, the linear system is solved twice at each time step. First with the 
matrices for the evaporative flux without coalescence (i.e., A3 = Bj = Cj = 0), and second 
with the matrices for coalescence without evaporation (i.e., E^ — 0). The unknowns a, . . ., 
C3 are found by adding the two solutions. 

As discussed earlier, for the evaporative step the linear system can be decomposed into five 
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NxN systems that can be solved sequentially. Likewise, for the coalescence step a and b* 
can be found separately by solving a 2Nx2N system: 



Ai 


A2 




a 




Pa 


A3 


A4_ 




b* 







(47) 



and then each of the vectors Cj can be found separately: 

Cjh*, 



B,a 



(48) 



where (for coalescence) Dj = V is a Vandermonde matrix [31] formed from the volume 
abscissas Other choices of moments have also been found to be numerically stable. For 
example, another valid choice is 



k^{i-l)/3 ie {!,..., 2N} 

A;=(2i-l)/3 ie{l,...,N} 

k^{2i-l)/3 ie{l,...,N} 

A;=(2i-l)/3 ie{l,...,N} 



with I — m — p — 

with I — 1, m — p — 

with m = 1, I — p — 

with p — I, I — m — 0. 



(49) 



This choice can be found to give more accurate results for our test cases and still includes the 
surface area and the volume of the droplets, as well as their momentum. Thus, it will be used 
for the computations in Section 4. We should note that for a given value of A^, the simulation 
results found with the moment sets in Eqs. (45) and (49) were nearly identical. The choice 
between these two systems was thus made based on ease of solution of the linear system. 

In general, when moments involving the velocity are limited to first order, the matrices that 
must be inverted will be non-singular as long as the volume abscissas are distinct. The nu- 
merical treatment of the singularities associated with Eq. (47) has been discussed elsewhere 
[24]. The coalescence operator will normally force the Vn to remain distinct if they have dis- 
tinct velocities. However, if due to initial conditions two or more of the volume abscissas are 
equal, it suffices to perturbate the values of Vn enough to allow for the coefficient matrix in 
Eq. (47) to be invertible. We should also note that for cases dominated by coalescence (e.g., 
without evaporation) the volume abscissas grow rapidly, leading to matrices that are more 
and more ill-conditioned as the abscissas increase. Thus, even though they are strictly non- 
singular, such matrices lead to severe numerical difficulties. However ill-conditioning can be 
almost completely alleviated by using iterative improvements of the linear solver as described 
in Section 2.5 of Press et al. [31] after rescaling Eq. (30). The latter is done by defining pos- 
itive scaling factors and Ug, and dividing both sides of Eq. (30) by v^uI^"^~^p. Note that 
the abscissas and unknown source terms are rescaled in a consistent manner: v-n Vn/vg, 
u„ — > Un/us, a ^ a, b* ^ h*/vs, and c* — > Cj/{vsUs). The evaporative flux constraints 
(Eqs. (35), (38) and (40)) can be rescaled in a similar manner by introducing a positive scaling 
factor Wg for the weights: w„ Wn/wg. In this work, we use the following scaling factors: 
Vg = max„ Vn, Ug = max„ |u„| and Wg = J2n We find that using the scaled variables and 
at most three iterative improvements of the linear solver are enough to completely eliminate 
round-off error in the solution to the DQMOM linear system. Moreover, because round-off 
error leads to poor performance of the differential equation solver, the overall computational 
cost using the iterative improvements can be significantly reduced. 
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3 Nozzle Test Problem 



In order to validate the proposed DQMOM approach for strongly coalescing sprays, and to 
compare this method to both a reference Lagrangian solver solution as well as the solution ob- 
tained with the multi-fluid model, we need a well-suited test problem that is difficult enough 
to highlight the limitations of the methods under consideration. For that purpose, we have 
chosen for the gas phase a 2D axisymmetrical conical decelerating nozzle, designed in such 
a way that it admits, for one-way coupling spray dynamics a self-similar solution. After pre- 
senting the details of this configuration, we will provide the set of DQMOM equations to be 
solved in this framework. We have selected six representative test cases, combining coales- 
cence/no coalescence with evaporation/no evaporation, which are then presented. Next we 
give an overview of the Lagrangian solver that provides the reference solution for the various 
test cases. Because the problems under consideration can be difficult to solve numerically, we 
must be very careful as far as this reference solution in concerned and thus we provide the 
details of the Lagrangian numerical integration in the limit of one-way coupling with the gas 
phase. Finally, before discussing the results in Section 4, we review the fundamentals of the 
Eulerian multi-fluid model for the sake of self-consistency of the paper. 

3. 1 Definition of configuration 

The chosen configuration is stationary 2D axisymmetrical in space and ID in droplet size. 
It is described in detail, along with the Lagrangian solver, in [22]. Hence, only its essential 
characteristics are given here. 

A spray of pure heptane fuel is carried by a gaseous mixture of heptane and nitrogen into a 
conical diverging nozzle of axis (0 < z). At the entrance, z = zq, 99% of the mass of the 
fuel is in the liquid phase, whereas 1% is in the gaseous mixture. The temperature (400 K) as 
well as the composition of the gas mixture (mass fraction is 2.9% for heptane and 97.1% for 
nitrogen) is fixed during the entire calculation. The gas density is then 0.871733 mg.cm"^. 
The influence of the evaporation process on the gas characteristics is not taken into account in 
our one-way coupled calculation. It is clear that the evaporation process is going to change the 
composition of the gas phase and then of the evaporation itself. However, we do not attempt to 
achieve a fully coupled calculation, but only to compare two ways of evaluating the coupling 
of the dynamics, evaporation and coalescence of the droplets. It has to be emphasized that it 
is not restrictive in the framework of this study, which is focused on the numerical validation 
of Eulerian solvers for the liquid phase under conditions of strong coalescence. 

For the problem to be one-dimensional in space, conditions for straight trajectories are used 
and are compatible with the assumption of an incompressible gas flow. This leads to the 
following expression for the gaseous axial velocity Vz and the reduced radial velocity Vr/r: 

,, = V(.) = '^, = £/(.-) = m = iZM to,>,„ (50) 
z'' r z z^ 

where > is the coordinate of the nozzle entrance and the axial velocity ^(^o) at the 
entrance is fixed. The trajectories of the droplets are also assumed straight since their injection 
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Figure 1. Initial number density functions for droplet radius. Left: Monomodal distribution. Right: 
Bimodal distribution. 



velocity is co-linear to the one of the gas. This assumption is only valid when no coalescence 
occurs. However, even in the presence of coalescence, it is valid in the neighborhood of the 
centerline. 



Let us finally consider two droplet distribution functions. The first one, called monomodal, 
is composed of droplets with radii between and 35 /xm, with a mean radius of 12 /xm, 
a variance of 5 fim and a Sauter mean radius of 15.6 fim. It is represented in Fig. 1 and 
is typical of the experimental conditions reported in [23]. The droplets are constituted of 
liquid heptane, their initial velocity is the one of the gas, their initial temperature, fixed at the 
equilibrium temperature 325.4 K (corresponding to an infinite conductivity model), does not 
change along the trajectories. The second distribution is called bimodal since it involves only 
two groups of radii, respectively, 10 and 30 microns with equal mass density. This bimodal 
distribution function is typical of alumina particles in solid propergol rocket boosters [17], 
and is represented in Fig. 1. 



The initial injected mass density is taken as mo = 3.609 mg.cm"^ so that the volume fraction 
occupied by the liquid phase is 0.57%. Because of the deceleration of the gas flow in the 
conical nozzle, droplets are going to decelerate, however at a rate depending on their size and 
inertia. This will induce coalescence. The deceleration at the entrance of the nozzle is taken as 
a(zo) = —2V{zq)/zq\ it is chosen large enough so that the velocity difference developed by 
the various sizes of droplets is important. We have chosen rather large values, as well as strong 
deceleration, leading to extreme cases: V{zq) = 5 m.s~^, zq = 10 cm for the monomodal case 
and V{zq) = 5 m.s^^, Zq = 5 cm for the bimodal case. These values generate a very strong 
coupling between coalescence, evaporation and droplet dynamics. These severe conditions 
as well as the two types of size distributions make the test cases under consideration very 
efficient tools for the numerical evaluation of the two Eulerian models. 
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3.2 DQMOM model equations in nozzle configuration 



For the nozzle test case, Eqs. (9-1 1) reduce to a set of ordinary differential equations (ODEs) 
defined on the interval z G [zq, oo) for the variables Wn, Vn, Vniz) = Ur/r and Cn(^) = u/. 

'^WnVn + (WnU) = ^n, (51) 
n) ^ni 

(52) 

'^WnVn'nl + {WnVnVnD = Cm/r, (53) 
^WnVnTlnin + dz (w„U„^^) = C^„, (54) 

where Uz = i{z) and Ur — rr]{z) are the axial and radial components of the spray velocity, 
respectively. The corresponding fluid velocities are given in Eq. (50). The terms on the right- 
hand side of Eqs. (52-54) are given by 

K = K + WnRv{Vn), (55) 
Crn/r = C*„ + WnllnRviVn) + tt^n^n^^r (^n, rin)/r, (56) 
C^n = C*„ + WninRv{Vn) + WnVnFz{Vn, C„), (57) 

where the drag model is 

Fr{v,r))/r^a(^—j {U-ri), (58) 
F,{v,0 = al^-j {V-0 (59) 

with a — 1.566 x 10"^ m^.s~^ 

From the form of the governing equations, it is straightforward to show that if rjn = inl^ 
at z — zq, then this relation will hold for all z and the droplet trajectories are straight lines. 
The system of DQMOM model equations can thus be reduced to three nonlinear ODEs for 
= Wn{z/zof, Vn, and by eliminating Eq. (53): 

dz (WnCn) = On, (60) 
dz (ty>nCn) =K + W^RviVn) (61) 



and 



dz {w:vnC) = + w*JMvn) + «<^n j {V - Q. (62) 

The terms on the left-hand side represent changes in the weights and abscissas due to trans- 
port. The terms on the right-hand side represent, respectively, the changes due to coalescence. 
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evaporation, and drag. The coalescence terms are found by solving 



N N N 

n=l n=l n=l 

1 fzo-' ^ 



^ ^ n=l 0=1 



(63) 



Note the presence of the scaling factor {zo/zf in the coalescence rate. As discussed in Sec- 
tion 2.4, we will use moments given in Eq. (49) that decouple Eq. (63) into two smaller 
systems. 



3.3 Test cases 



For evaporation, we will consider three cases described below: (i) no evaporation (i?^ = 0), 
{ii) linear evaporation {R^ oc v), and (///) non-linear evaporation {Ry oc v^/^). For each 
case, we will consider two sub-cases: without coalescence (.Bcoai — 0) and with coalescence 
(^-coai = !)• The two evaporation laws correspond to the two cases described in Section 2.3, 
for which droplets disappear either in infinite time {ii), thus leading to a evaporative flux at 
zero size, or in finite time {Hi) for which the evaporative flux depends on the structure of the 
number density function in size phase space. 



No evaporation 

For the special case of no evaporation and no drag, the right-hand sides of Eqs. (52-54) are 
null. This special case has an analytical solution with , and ^„ constant. In the opposite 

limit of no evaporation and infinite drag, ^n — ^ and w* oc (^/^o)^• 
For non-evaporating droplets, — 0. In the absence of coalescence, a„ = 6* = c*„ — 0. 
The DQMOM model reduces to v„ and constant, and 

en9.en = «(^ — j (y-Cn). (64) 

This result is consistent with our earlier remark concerning the cases of zero and infinite drag. 
Finally, we should note that even with coalescence the momentum is conserved {k = m = l) 
so that c*zn — 0- Thus, we can expect to be approximately constant for all values 
of drag. For this case we expect excellent agreement between DQMOM and the Lagrangian 
solver in the absence of coalescence since the corresponding transport equations are identical 
(i.e., each DQMOM abscissa behaves like a Lagrangian particle). On the other hand, with co- 
alescence the droplets grow very large and we expect differences due to how the coalescence 
term is treated in each method. This test case will, however, be very difficult for the multi- 
fluid model, since it was especially designed to tackle the problem of evaporation. In the 
presence of strong growth of droplet size, the number of sections that must be used in order 
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to reproduce the physics with the multi-fluid model will also dramatically increase. Conse- 
quently, this test case will allow us to both test the capability of the DQMOM to capture the 
coupling of dynamics and coalescence at low cost in comparison to the Lagrangian solution, 
and to see if the multi-fluid model can provide good results, even if it is not competitive in 
terms of computational efficiency. 

Linear evaporation 

For evaporating droplets with linear evaporation, we take 

Rv{Vn) = -EyVn, (65) 

with Ey — 7.1262 s~^ for the monomodal case and = 14.2524 s~^ for the bimodal case. 
For this case, the evaporative flux ^ is zero. The coalescence terms are again found by solving 
Eq. (63). In the absence of coalescence, we have a„ = 6* = c*„ = and the DQMOM model 
reduces to w*^„ constant, Eq. (64), and 

^ndzVn = Rv{Vn)- (66) 

Thus the volume Vn and velocity are coupled through evaporation and drag, but are inde- 
pendent of w* in the absence of coalescence. For this case we again expect excellent agree- 
ment between DQMOM and the Lagrangian solver in the absence of coalescence since the 
corresponding transport equations are identical. On the other hand, with coalescence there is 
a competition between growth and evaporation leading to smaller droplets than in the non- 
evaporating case. This is a very interesting test case, since it will allow us to compare both 
methods in an evaporative configuration, but without getting into the difficulty of modeling 
the droplet disappearance with the DQMOM approach. 

Non-linear evaporation 

With non-linear evaporation we will use 



with Eg = 1.99 X 10~^ m^/s. For this case the evaporative flux ijj will generally be non-zero, 
and is found using the method described in Section 2.3 with w;* in place of Wn- However, 
we will also compare predictions for the bimodal initial distribution found by setting ^ = 0. 
As for the previous cases, we will investigate the effect of the flux model with and without 
coalescence. From a practical standpoint, the behavior of DQMOM with non-linear evapora- 
tion is of great interest and it is a configuration with which the comparison of both Eulerian 
models will be of practical relevance. 

3.4 Reference Lagrangian solution 

Euler-Lagrange numerical methods are commonly used for the calculation of polydisperse 
sprays in various application fields (see for example [30,18,27,28,33,10] and the references 




(67) 
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therein). In this kind of approach, the gas phase is generally computed using a deterministic 
Eulerian solver, while the dispersed phase is treated in a Lagrangian way. The influence of the 
droplets on the gas flow is taken into account by the presence of source terms in the system 
of gas conservation equations. Two Lagrangian methods can be used as far as the dispersed 
phase is concerned depending on the level at which the physical processes are modeled. The 
first one is a Discrete Particle Simulation, which describes the evolution of numerical parti- 
cles, each one representing one or several droplets. The physical processes such as transport, 
evaporation, collisions are then described by Liouville equations and the Eulerian fields usu- 
ally recovered through ensemble averages. However, in the present study, we have preferred 
the Williams governing equation and thus a statistical description of the coalescence process. 
We then coherently use a Direct Simulation Monte-Carlo method (DSMC), the second kind 
of approach. It can be seen as the uncoupling, over a small time step, of the droplet transport 
in phase space (dynamics and evaporation), described by a particle method, and the collisions 
described by a Monte-Carlo method. 

A complete exposition on the derivation and implementation of this method is outside the 

scope of this paper. We refer the reader, for example, to [2,18,17] for more details. Here, for 
the sake of completeness, we present only the main features of the numerical method that we 
used in order to provide a "reference numerical solution" for the test cases. 



Lagrangian solver 

The Lagrangian solver can be roughly interpreted as a stochastic representation of the kinetic 
equation (1). In other words, in the limit of a sufficiently large number of stochastic parti- 
cles and a sufficiently fine computational grid (at least in the case of one-way coupling), the 
statistical estimates for the moments found from the particles should converge to those com- 
puted from the Eq. (1). In the Lagrangian solver, at each time step k, the droplet distribution 
function f{t'') is approximated by a finite weighted sum of Dirac masses, f{t'^), which reads 

i=l 

Each weighted Dirac mass is generally called a "parcel" and can be physically interpreted as 
an aggregated number of droplets (the weight nf), located around the same point, zf, with 
about the same velocity, and about the same volume, v^. iV^ denotes the total number of 
parcels in the computational domain at time t'^. In all our calculations, the weights were 
chosen in such a way that each parcel represents the same volume of liquid (nf = Const). 

Each time step of the particle method is divided in two stages. The first is devoted to dis- 
cretization of the left-hand side of the kinetic equation (1), modeling the motion and evapo- 
ration of the droplets. In our code, the new position, velocity and volume of each parcel are 
calculated according to the following numerical scheme: 



1 exp(- Ai/rf ) + V{z^) (l - exp(-AVTf )) (69) 
dv/R^{v) = At (70) 

k 
i 

z^+^ = + At u'l+^ = + At V{z^,) + At (uf - Viz'l)) exp(-AVrf ) (71) 
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where V denotes the axial gas velocity, (u'^) corresponds to the axial coordinate of the 
position (velocity) of the parcel i at time t'^, Rv is the evaporation rate (independent of t and 
X since the gas composition and temperature are assumed constant in the domain). Eq. (70) is 
resolved analytically and depends on the chosen evaporation model. For linear evaporation, 
it can be written as 

v^+^ ^v^eM-Ev^t) (72) 
and for non-linear evaporation, it is written as 

sf+i = - EsAt (73) 
where sf is the parcel surface area. The parcel relaxation time is defined as 

r' = ^ (74) 

with rf being the parcel radius, p the liquid density and /ig the gas viscosity. 

In Eqs. (69-71), the parcel radial coordinate is not calculated because the trajectories of the 
parcels are straight lines. Besides, as mentioned above, the influence of the droplets on the 
gas flow is not taken into account. Hence, Eq. (50) is used to calculate the gas velocity, V{z^), 
at the parcel location. 

The second stage of a time step is devoted to the discretization of the collision operator. 
Several Monte-Carlo algorithms have been proposed in the literature for the treatment of 
droplet collisions [30,18,33,17,34]. They are all inspired by the methods used in molecular 
gas dynamics [4] and, in particular, they suppose that the computational domain is divided 
into cells, or control volumes, which are small enough to consider that, within them, the 
droplet distribution function is almost uniform. 

The algorithm used in our reference Lagrangian solver is close to the one proposed by 
O'Rourke. It consists of the following 3 steps (see also [18] for more details): 

1. For each computational cell Cj, containing Nj parcels, we choose randomly, with a uni- 
form distribution law, Nj/2 pairs of parcels, {Nj — l)/2 if Nj is odd. 

2. For each pair p, let pi and p2 denote the two corresponding parcels with the convention 
rii > n2, where rii and n2 denote the parcel numerical weights. Then for each pair p of the 
cell Cj, we choose randomly an integer Up, according to the Poisson distribution law: 

-P(i^) = ^exp(-Ai2), 



with 



rH{Nj-l)At 2| I 

Al2 = TT — — + R2) Ml - M2 

vol(Cj) 



with vol(C,/) being the volume of cell Cj, which is proportional to {zj/zq)"^ for the nozzle 
test case, and Ri, R2 being the radii of the parcels pi,p2 - The coefficient A12 represents the 
mean number of collision, during (Nj — 1) time steps, between a given droplet of parcel p2 
and any droplet of parcel pi. Note that a given pair of parcels is chosen, on average, every 
{Nj — 1) time steps. 
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If Up — 0, no collisions occur during this time step between the parcels pi and p2. Oth- 
erwise, if Vp > 0, parcel pi undergoes Up coalescence with parcel p2 and the outcome 
of a collision is treated as follows. First the weight rii of the parcel pi is replaced by 
n[ = Til — ypn2 and its other characteristics are left unchanged. If n'^ < 0, parcel pi is 
removed from the calculation. Secondly, the velocity U2 and the volume V2 of parcel p2 are 
replaced by 



V2 = V2 + I^pVl, U'2- ■ , 

V2 + fpVi 

and its weight, 712, is left unchanged. 

Let us mention that, for each time step and each control volume Cj, the computational cost 
of this algorithm scales like 0{Nj). This is a great advantage compares to the O'Rourke 
method, which scales like O(A^j). Another algorithm, with the same scaling features, has 
been introduced by Schmidt and Rutland in [34]. 

To obtain good accuracy, the time step, Ai, must be chosen small enough to ensure that the 
number of collisions between two given parcels, pi and p2, is such that for almost every time: 
Vpn2 < rii. The average value of Up being A12, this constraint is equivalent to the condition 



, V2U2 + VpViUi 



n2NjAt_^^ , „ 2 



vol(Cj) 



7r{Ri + R2r\ui - U2\ ^ 1. (75) 



For the nozzle test case described above, this constraint reveals to be less restrictive than the 
"CFL" condition 

Vi = l,....7V, (76) 

with Az being the mesh size. This condition is necessary to compute accurately the droplet 
movement and in particular to avoid that a parcel goes through several control volumes during 
the same time step. This is essential in order to have a good representation of the droplet 
distribution function in each mesh cell. 



Reference solution 

The Lagrangian solver just described is used to provide reference solutions in stationary cases 
with and without coalescence. In order to obtain a converged solution, particular attention 
must be devoted to the choice of the number of parcels, the size of the cells, and the time 
step. 

For cases without coalescence, the computational cells are only used to have spatially av- 
eraged quantities to compare with Eulerian results. Moreover, the stationary aspect of the 
problem allows averaging in time in order to obtain smooth solutions. For these reasons, 
the conditions on the number of parcels and on the size of the computational cells are not 
very restrictive in the absence of coalescence. The time step is only limited by the CFL-like 
condition (76) needed for the convergence, with a low value. This last condition is the most 
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restrictive since the scheme used for the transport of the particles is first order. For our test 
cases, the time step must be 10"^ s or smaller. 



distribution 


evaporation 


No. of parcels 


No. of parcels inj./s 


monomodal 


linear 


41,560 


100,000 


monomodal 


nonlinear 


20,440 


1,000,000 


bimodal 


nonlinear 


6,320 


200,000 



Table 1 

Number of parcels for the Lagrangian simulations for the cases without coalescence. 



distribution 


evaporation 


No. of parcels 


No. of parcels inj./s 


min. No. of 
parcels /cell 


monomodal 


linear 


160,000 


200,000 


40 


bimodal 


linear 


126,000 


560,000 


50 


monomodal 


nonlinear 


35,000 


1,300,000 


260 


monomodal 


no 


44,200 


300,000 


65 



Table 2 

Number of parcels for the Lagrangian simulations for the cases with coalescence. 

For cases with coalescence, there are additional restrictions. First, the algorithm used for co- 
alescence assumes that the droplet distribution function of the spray is nearly uniform over 
each computational cell. However, in the region with high gradients of the gas velocity, that 
is to say at the entrance of the nozzle, this distribution can change quickly and the size of 
the cells must be small enough to avoid numerical errors. Moreover, in order to properly de- 
scribe the coalescence phenomenon in each cell with the stochastic algorithm, a sufficient 
number of parcels must be present in each cell, typically on the order of 50, with a mini- 
mum of 20 [1]. The smaller are the cells, the larger must be the number of parcels in the 
computational domain. The required size of the cells is evaluated for the case where the size 
distribution changes the most rapidly (the case without evaporation). We then employ a non- 
uniform space discretization with 130 cells, with smaller cells near the entrance of the nozzle 
defined using a uniform discretization for the variable z^^^^ . The number of parcels injected 
per second is given in Tables 1 and 2. 

3.5 Eulerian multi-fluid solver 

Eulerian multi-fluid methods were developed as an alternative to Lagrangian methods for 
the simulation of poly disperse evaporating sprays. A complete derivation of such methods 
from the kinetic model is performed in [21] for dilute sprays and in [22] for sprays with 
coalescence. The principle of the method is quite different than the one used in DQMOM. 
Indeed, it can be considered as a finite-volume discretization in the droplet size phase space 
for moments of order and 1 of the velocity distribution conditioned on size. 
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In laminar flows, it can be proven rigorously that it is sufficient to work with only these two 
moments as long as the velocity distribution conditioned on droplet size is mono-kinetic [7,8] 
(i.e. all droplets with the same volume have identical velocities so that the size-conditioned 
velocity variance is null.) By construction, the nozzle test problem will be mono-kinetic for 
non-coalescing droplets. However, with coalescence the size-conditioned velocity variance 
can be nonzero. Comparisons between the Lagrangian and Eulerian results in the presence of 
coalescence will therefore allow us to quantify the magnitude of the error caused by invoking 
the mono-kinetic assumption in the Eulerian models. (Recall that the choice of moments used 
in the DQMOM linear system is equivalent to the mono-kinetic assumption in the multi-fluid 
model.) In this section, we provide only the main points of the derivation of the multi-fluid 
model, as well as the underlying assumptions that are implied, and the resulting system of 
equations that will be solved. 

The first step consists of writing equations for the two moments in velocity. This leads to the 
closed semi-kinetic model if the following assumption is made concerning the structure of /: 
f{v, u; X, t) — n{v; x, t)S{u — u{v;x, t)). In other words, the droplet velocity conditioned 
on the size is assumed to be Dirac delta function. In the case of a coalescing spray, the 
compatibility of such a condition with droplet coalescence is far from obvious; however, the 
semi-kinetic system of conservation equations can be obtained by using an asymptotic limit 
as presented in [22]. 

The second step consists of discretizing n{v) in sections [v^^^^\v^^'') and in integrating the 
semi-kinetic model over each section. This leads to a multi-fluid model (by using a presumed 
distribution k,^^^ (v) in each section), thereby yielding a conservation equation on the moment 
associated with the mass density 

n{v;:x.,t) — m^^\x.,t)K^^\v) where / pvK^\v)dv — 1. 

In addition, only the averaged velocity is considered in each section, i.e. u{v; x, t) = u^-^^ (x, t) 
if f < V < v^^\ The resulting system can be found in [22]. It can be rewritten and sim- 
plified in the stationary, self-similar, 2D axisymmetrical configuration we are considering. 

The resulting set of equations is 

(i) 



- (E?) + Ei'^)m^^W^ + Ep+')m(^'+i)«,(^'+^) + m^'^F^^^ + C^), (78) 

where rt^^^^ is the axial velocity, which only depends on z, and ruj^^^ /z is the radial ve- 
locity, since the trajectories are straight lines. Moreover, E^^ and E^^ are the "classical" 
pre-calculated vaporization coefficients [13,21]: 

E? = -pv^^-^'^ Rviv^'~^^) K^'\v^'-^h and E^'^ = - / pR^{v) k^'\v) dv, 
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) is the axial component of the "classical" pre-calculated drag force [13,21]: 



F^O) = / pvF,{v,u,^^'>) K^^\v)dv where v^^^ ^ 

Ji;0-1) 



ji%vy^KU){v)dv^ 



The source terms associated with coalescence phenomenon in the mass and momentum equa- 
tion, respectively, of the jth section read 

N I'-J'i 
k=l 1=1 



AT 
k=l 



lU) 



i=l 



where Vjk — \uz^' — uj^ ' \ and the collision integrals Qjk, Q% and Q*^ do not depend on z. 
The disappearance integrals Qjk are evaluated on rectangular domains Ljk — [v''^~^\ v'^^^] x 
,y(fe)j^ whereas the appearance integrals, and are evaluated on the diagonal 
strips = {(f*, V*), v''^^^^ < if + V* < v'^^'>} / ^^^=1 ^kk, which are symmetric strips with 
respect to the axis v"^ = v*. These strips Dj* are divided into domains, denoted by Xji and 
the symmetric one, XjJ^, where the velocity of the partners is constant. The domains Xji 
and Xjf^ are the intersection of Dj* with L^i, k > I, and L^i, k < I, respectively; their 
index is denoted i e [1, Z*^^^] and we define two pointers that indicate the collision partners 
for coalescence, at fixed i: Oj^ — k and o*^ — l. 

The coefficients used in the model, either for the vaporization process or the drag force Ei\ 
Ei^^ and F^^\ j = [1,N] in Eqs. (77-78), or for the coalescence: Qj^, J = [l,N],k = 
[1, N],k j, Q]i, Q*i, j = [2,N],i = [1, /(j)] in Eqs. (79-80) can be pre-evaluated from 
the choice of k,'^^ in each section. The algorithms for the evaluation of this coefficients are 
provided in [22] . The distribution function is chosen constant as a function of the radius in the 
sections 1 to and exponentially decreasing as a function of the surface in the last section, 
as done in [22]. 

Because only the one-way coupled equations are solved and since the structure of the gas 
velocity field is prescribed and stationary, we only have to solve the ID ordinary differential 
Eqs. (77,78) for each section. The problem is then reduced to the integration of a stiff initial 
value problem from the inlet where the droplets are injected until the point where 99.9% of the 
mass has evaporated. The integration is performed using LSODE for stiff ordinary differential 
equations from the ODEPACK library [16]. It is based on BDF methods [14] (Backward 
Differentiation Formulae) where the space step is evaluated at each iteration, given relative 
and absolute error tolerances [16]. The relative tolerance, for the solutions presented in the 
following, are taken to be 10"^, and the absolute tolerance are related to the initial amount 
of mass in the various sections, since it can vary of several orders of magnitude. Repeated 
calculations with smaller tolerances have proved to provide essentially the same solutions. 
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4 Results and Discussion 



Simulations for the cases presented in the previous section were carried out with the La- 
grangian method, the multi-fluid model, and DQMOM. Except for the cases without evap- 
oration for which the multi-fluid method is not well suited (it requires a large number of 
sections and is only presented for comparison purposes), the Eulerian methods were solved 
using a initial-value solver for ODEs and required very short computational times (i.e., CPU 
sees) on a desktop computer. It is interesting to note that in the case without evaporation, the 
small computational cost still holds for the DQMOM approach. 

In contrast, the time- and ^-dependent Lagrangian simulations required several CPU hrs for 
each case. Because the DQMOM and multi-fluid results do not depend on time, it is not ap- 
propriate to compare the computing times directly. Nevertheless, it will generally be the case 
that using Eulerian methods will result in a substantial reduction in the computing time for 
solving the spray equation. Such a statement was studied in details in [22] for unsteady calcu- 
lations and the conclusions drawn from that paper are applicable to the two Eulerian methods 
presented here. Thus, the principal open question is whether or not the DQMOM results are 
of comparable accuracy to the multi-fluid model and to the more costly Lagrangian simu- 
lations. We will compare predictions for selected statistics from the three solution methods 
in order to answer this question. For the monomodal distribution and DQMOM resolution, 
several initial conditions will be used in the following an are presented in Table 3. 



Monomodal distribution 




Vol. moments, N=A 


Rad. moments, N=4 


Rad. moments, N=6 


Rad. moments, N=S 


n 




rn 




rn 


Wn/No 


rn 


Wn/No 


rn 


1 


0.7323 


9.9955 


0.1845 


4.4079 


8.5573E-2 


3.3423 


4.6445E-2 


2.8465 


2 


0.2545 


18.5282 


0.5397 


11.0409 


0.2779 


7.5262 


0.1488 


5.5373 


3 


1.288E-2 


27.5630 


0.2635 


18.2840 


5.5339E-2 


12.9743 


0.3089 


9.6916 


4 


2.279E-4 


36.0142 


1.212E-2 


28.3910 


4.9778E-3 


18.8823 


0.3438 


14.2697 


5 










3.1137E-4 


26.3693 


0.12931 


19.2986 


6 










1.6671E-5 


34.7171 


2.0905E-2 


25.2866 


7 














1.6982E-3 


31.5808 


8 














6.5627E-5 


37.5149 



Table 3 



Initial conditions for weights and abscissas found using QMOM. 



The representative moments used to compare the three solution methods are the number den- 
sity mo, the mass density mi, the average axial velocity difference between droplets and the 
gas phase Ud, and the Sauter mean radius rs2- They are defined by 

mo — //(v, u)dvdu, mi— pvf{v, u) dv du, 
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1/3 / ''./(''■ ^) 

/ v'^l^f{v, u) dvdu' 



With the DQMOM approach, these quantities are written 



N N 



n=l n=l 



N 



l/3^n=l^»W 

v^AT 2, 
E„=l 



= m ^nVniin - V)/mi, 



r32 = (47r/3) 



2/3 ■ 



'n 



n=l 

And with the multi-fluid method, they are 




N 

nil = 

i=i 



= ^m^-'''(M2*--^-* — V)/mi, 



N 



rz2 = (47r/3) 



Note that in practical applications, the mass density is a key quantity because it represents the 
total mass of liquid contained in the droplets. In the nozzle test case, the rate of coalescence is 
strongly dependent on the velocity difference between droplets, which we find to be strongly 
correlated with the average axial velocity difference. Indeed, if is not accurately captured, 
then we find that the predictions for all moments will degrade accordingly. In addition to the 
moments, we will also compare the mean droplet velocity conditioned on the radius {uz\r) at 
selected downstream locations, as well as the mass distribution function (pvf). For the DQ- 
MOM, the scaled weights will be used to represent the mass distribution function. Obviously, 
since the sum of the weights equals the area under the mass distribution function, the absolute 
value of the heights of the scaled weights is arbitrary. Nevertheless, the relative heights and 
the locations provide insight into how well the quadrature points represent the distribution 
function. 

We should note that for the monomodal cases without coalescence, the results with no evap- 
oration were essentially identical for all three solution methods. The results presented below 
for the monomodal case with linear evaporation are representative of the quality of the pre- 
dictions for all cases without coalescence and no evaporation. Likewise, for the bimodal case 
without coalescence and with linear or no evaporation, DQMOM and the Lagrangian method 
were essentially identical. The multi-fluid method also yielded very good results for these 
cases if the number of sections was chosen large enough to mitigate the numerical diffusion 
in the size phase space associated to the description of evaporation that leads to broadening of 
the peaks. Nonetheless, because none of these cases revealed any unanticipated problems for 
any of the simulation methods, we will not discuss them further. Instead, we will primarily 
focus on cases that present particular challenges to one or more of the solution methods. 

4.1 Monomodal case: linear evaporation without coalescence 

We begin with a representative case where all three solution methods yield essentially iden- 
tical results for all statistics. As noted in the discussion of the methods, for linear evap- 
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radius (microm) radius (microm) 



Figure 2. Monomodal case with linear evaporation {R^ 7.126?;). Top left: mass density. Top right: 
velocity difference. Bottom left: mass distribution function {z = 16 cm). Bottom right: conditional 
velocity {z = 16 cm). 

oration without coalescence the DQMOM equation for each node is the same as the La- 
grangian model. Thus, the only difference between the two solution methods is that the La- 
grangian method uses many more particles to represent the spray than the DQMOM method 
(A^ = 4). For the monomodal distribution, the multi-fluid model does not require many sec- 
tions (N = 10) to accurately capture cases with linear evaporation without coalescence. The 
simulation results for the three methods are shown in Fig. 2. It can be observed that the mass 
density and velocity difference predicted by the three methods are nearly identical. From the 
plot of the mass distribution function at 5; = 16 cm, we can see that the multi-fluid model 
with ten sections does a good job of capturing the Lagrangian mass distribution function. 
Likewise, the DQMOM weights and abscissas follow the general shape of the Lagrangian 
mass distribution function. Finally, for the conditional velocity {uz\r) we see that all three 
methods produce the same curve. We should note that for cases without coalescence the La- 
grangian simulations predict essentially no velocity dispersion about the conditional value. 
In other words, conditional velocity fluctuations defined by u'{r) = {{uz — {uz\r))'^\ry^'^ are 
null. This is exactly one of the necessary conditions evoked when deriving the multi-fluid 
model, which would explain why its predictions for this case are in excellent agreement with 
the Lagrangian method. 

4.2 Monomodal case: nonlinear evaporation without coalescence 

Cases with nonlinear evaporation result in a loss of droplets in finite time, which translates 
into a nonzero flux ipit) in DQMOM. For the monomodal case without coalescence, we 
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Figure 3. Monomodal distribution function with optimal sections. 



expect the flux term to be a smooth function of t, and thus it cannot be neglected. In the 
multi-fluid model, the flux is computed directly from the shape of the first section (i.e., the 
section near the origin) and does not yield any difficulty. 

In our multi-fluid simulations, we use the "optimal" choice of sections with = 12 shown 
in Fig. 3 [23]. Obviously, a finer discretization (larger N) could be used in the multi-fluid 
model to attain closer agreement with the Lagrangian method, but this would increase the 
computational cost. Note that the first section is represented by a constant slope, which cor- 
responds to a constant flux level at each time step. For the DQMOM, we use = 4 and 
the evaporative flux is computed using the ratio constraints introduced in Section 2.3. It can 
be noticed that the increase of N do not imply an increase of the number of conserved mo- 
ments during the evaporation step since the number of ratio constraints is also increasing in 
the same way. The value of N is then conditioned by the capacity of the method to follow the 
dynamics of droplets of different sizes. Representative results for the three solution methods 
are shown in Fig. 4. In general, all three methods produce very similar predictions. From the 
number density, we can observe that DQMOM with the ratio constraints does a good job of 
predicting the loss of droplets due to evaporation. Likewise, the mass densities found from all 
methods are very close. We should note that for z > 20 cm the number of remaining droplets 
is very small and the statistics computed from the Lagrangian method are subject to statistical 
errors. Comparing the Sauter mean radii predicted by the three methods, we can observe that 
the agreement is generally satisfactory up to 2; = 20 cm. The DQMOM shows the largest 
deviation from the Lagrangian method at 2; = 20 cm due to errors in the flux model, but 
the agreement is still acceptable. The differences in the Sauter mean radius are reflected in 
the predictions of the velocity difference. In general, droplets with a larger radius will have 
a higher velocity difference. Thus, we see that initially the Sauter mean radius predicted by 
DQMOM is larger than that from the Lagrangian method, resulting in a slightly higher ve- 
locity difference at 2; = 12 cm. Later on (2; > 15 cm) this trend is reversed. Finally, we can 
note that neglecting the flux in DQMOM yields poor predictions of number density since we 
can observe the artificial jumps in the number density related to the singular fluxes associated 
to one abscissa crossing the zero size limit, as well as the oscillating dynamics of the Sauter 
mean radius for this case. 
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Figure 4. Monomodal case with nonlinear evaporation. Top left: number density. Top right: mass 
density. Bottom left: Sauter mean radius. Bottom right: velocity difference. 



4.3 Bimodal case: nonlinear evaporation without coalescence 



By changing from the monomodal to the bimodal distribution, we change the nature of the 
initial distribution function and thus the nature of the numerical difficulties. For the multi-fluid 
model, the bimodal case is difficult because a relatively large number of sections (A^ = 30) 
is needed to capture the two peaks with acceptable numerical diffusion. The use of a second- 
order method developed in [20] would reduce this number to around 10; however, it would 
still be difficult to describe Dirac delta function by a finite-volume approximation. On the 
other hand, this case is "optimal" for DQMOM because only two {N = 2) abscissas are 
required (one for each peak) and the flux is null, expect when a peak passes the origin. In 
Fig. 5 results from the three simulation methods are shown and it is clear that DQMOM 
performs extremely well for this case by setting = 0. For example, the number density 
function shows step changes at z = 7.2 cm and 13.8 cm (i.e., when a peak passes the origin), 
and DQMOM exactly reproduces this behavior. With N = 30, the multi-fluid model does 
a good job of predicting the mass density. However, from the plots of number density and 
Sauter mean radius, we can observe the negative effects of numerical diffusion, which tends 
to smooth out the peaks in the distribution (the method has been shown to be first order 
in the droplet size discretization step in [20], where some higher-order methods have been 
proposed). Nevertheless, all three methods yield reasonable predictions for all of the cases 
without coalescence. We should note, however, that for more complicated initial distributions 
(e.g., delta functions combined with smooth functions) specifying the evaporative flux in 
DQMOM may be problematic. 
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Figure 5. Bimodal case with nonlinear evaporation. Top left: number density. Top right: mass density. 
Bottom left: Sauter mean radius. Bottom right: velocity difference. 



4.4 Monomodal case: linear evaporation with coalescence 



We now turn to the more difficult cases that include coalescence. As mentioned earlier, the 
coalescence of droplets with different volumes (and velocities) will lead to velocity disper- 
sion (n'(r) > 0). Physically, this implies that two droplets with the same volume will have a 
nonzero probability of colliding (due to the difference in velocity). Thus, the rate of coales- 
cence when u'{r) > will be larger than when u'{r) = 0. Numerical approximations (such 
as the multi-fluid model) that assume n'(r) = should therefore predict smaller droplets than 
the Lagrangian method. In Fig. 6 we present results for the three methods for linear evapora- 
tion (-ip = 0) with coalescence. ^From the velocity difference we can observe that coalescence 
leads to a slower relaxation to the gas velocity due to formation of larger droplets than without 
coalescence. Note that in general all three methods predict similar results for the velocity dif- 
ference. However, due the differences in the predictions of the shape of the mass distribution 
function, the multi-fluid model predicts slightly slower relaxation and the DQMOM slightly 
faster than is found with the Lagrangian method. Comparing with Fig. 1, we can observe that 
coalescence leads to much larger droplets than are present in the initial distribution function. 
In general, the multi-fluid model predicts a slightly larger number of droplets above 80 /im 
than the Lagrangian method. Nevertheless, the predictions are in reasonably good agreement. 
The predictions for the conditional velocity {uz\r) are also good. Finally, note that we used 
= 6 with DQMOM, the reason for which will be discussed for a more difficult case in 
Section 4.7. 
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radius (microm) radius (microm) 

Figure 6. Monomodal case with linear evaporation {Ry ^ 7.126?;) and coalescence. Top left: mass 
density. Top right: velocity difference. Bottom left: mass distribution function (z = 22 cm). Bottom 
right: conditional velocity (z = 22 cm). 

4.5 Bimodal case: linear evaporation with coalescence 

We now consider a more difficult case where the initial distribution is bimodal. As discussed 
previously, the peaks in the distribution are difficult to resolve accurately in the multi-fluid 
model with a limited number of sections . When combined with coalescence, this has impor- 
tant consequences because numerical diffusion can lead to spurious coalescence of droplets 
with slightly different volumes (and hence velocities) as observed in [22]. For example, with 
the bimodal distribution with droplets of radii 10 and 30 /im, coalescence cannot produce 
droplets below 30 jim. However, spurious coalescence between droplets of radii near 10 /im 
leads to droplets in the range below 30 /im. We overcome this difficultly by using a large 
number of sections (A^ = 500) in the multi-fluid model. This number could also be reduced 
by using a second-order method for the evaporation such as the one of [20] but this is not 
the point we want to make with this configuration. Note that the same problem arises in the 
Lagrangian method when the spatial cell size /S.z is too large. While DQMOM does not suffer 
from spurious coalescence, the bimodal case is still difficult because the initially two-peak 
distribution will quickly form multiple peaks due to pair- wise collisions. With = 6 in DQ- 
MOM, it is at best possible to represent six peaks. Results for the three methods are shown in 
Fig. 7 where it can be seen that the mass density and the velocity difference are reasonably 
well predicted by the multi-fluid model and DQMOM. From the mass distribution function 
at z = 11 cm, the multi-peak structure due to coalescence is quite evident, as is the slight 
numerical diffusion in the multi-fluid model (even with A^ = 500, but this is expected since 
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Figure 7. Bimodal case with linear evaporation {Ru 14.252t)) and coalescence. Top left: mass 
density. Top right: velocity difference. Bottom left: mass distribution function {z = 11 cm). Bottom 
right: conditional velocity (z = 11 cm). 

this is a first-order method). Note that DQMOM with = 6 has two abscissas at points cor- 
responding to the major peaks (10 and 33 fim), and the remaining abscissas at points without 
major peaks. Comparisons of the conditional velocity predicted by the three methods are also 
quite favorable for this difficult case. 



4.6 Monomodal case: nonlinear evaporation with coalescence 



We will now consider the more physically relevant case of nonlinear evaporation. As dis- 
cussed earlier, the evaporative flux for this case is nonzero, so we will need to model it in 
DQMOM. Here, we consider two models for '0: (a) ratio constraints and (b) if) = {i. Because 
the initial distribution is monomodal, we might expect that using ratio constraints is always 
a better choice. On the other hand, if coalescence is much faster than evaporation, it might 
happen that droplets grow faster than they disappear so that the evaporative flux is closer to 
zero. For the multi-fluid model, we use = 15 sections. Results for the three methods are 
shown in Fig. 8. The number density illustrates the effect of the choice of if) in DQMOM. 
With = 0, the number density changes discontinuously whenever an abscissa passes the 
origin. However, DQMOM with ratio constraints yields predictions very similar to the other 
two methods. Likewise, the mass density is predicted to be very similar for all three meth- 
ods; however, using zero flux with DQMOM is slightly worse. The predictions for the Sauter 
mean radius show opposing trends. In general, the multi-fluid model overpredicts the mean 
radius (i.e. predicts too much coalescence), while DQMOM underpredicts it. As before, for 
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Figure 8. Monomodal case with nonlinear evaporation and coalescence. Top left: number density. Top 
riglit: mass density. Bottom left: Sauter mean radius. Bottom right: velocity difference. 

the DQMOM predictions, the results with the ratio constraints are best. The predictions for 
the velocity difference follows the same trend. As discussed in the next example, the differ- 
ences observed between the Lagrangian method and the two Eulerian methods is likely due 
to the latter's inability to capture velocity dispersion. Moreover, we have used = 6 with 
DQMOM since, as shown in Section 4.7, it is adequate to describe coalescence phenomenon 
for this particular set of moments. 



4. 7 Monomodal case: coalescence with no evaporation 



In order to highlight the role of coalescence on determining the evolution of the number den- 
sity function, we now consider a case with no evaporation. For this case, droplets will grow 
continuously due to coalescence, and velocity dispersion will enhance the collision rate and 
lead to even larger droplets. Because the multi-fluid model uses fixed sections, it is necessary 
to fix the maximum radius at 200 /xm with N = 500 in order to capture the largest droplets 
at 2 = 30 cm. In contrast, the abscissas in DQMOM move in phase space to accommodate 
growth. Nevertheless, we can anticipate that the number of abscissas will affect the accuracy 
of the DQMOM predictions. In Fig. 10 it can be observed that when the number of moments 
increases, the accuracy of the DQMOM solution increases, from something almost ignoring 
the coalescence phenomenon with = 2 to a saturation of the accuracy for > 8. Indeed, 
the accuracy of the DQMOM for the description of the coalescence is related to the accuracy 
of the approximation of the coalescence operator by the quadrature formula (24) and (26), 
which increases with N. Since the results are quite good and at a low cost (and the linear 
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radius (microm) radius (microm) 

Figure 9. Monomodal case with coalescence but no evaporation. Top left: mass density. Top right: 
velocity difference. Bottom left: mass distribution {z = 22 cm). Bottom right: conditional velocity 
{z = 22 cm). 



system is reasonably well conditioned), we will use = 6 for comparisons with the other 
two methods. 

As mentioned earlier, without coalescence or evaporation all three methods predict essen- 
tially identical results. In Fig. 9 the results for the pure coalescence case are shown. Notice 
that the mass density does not decrease to zero because there is no evaporation; however, it 
does change due to transport. From the velocity difference, we can see that the multi-fluid 
model and DQMOM overpredict the relaxation rate. As discussed previously, this is due to 
both methods underpredicting the mean droplet size. From the mass distribution functions at 
z = 22 cm, we can observe that the Lagrangian method has more droplets with radii above 
80 /im than the multi-fluid model, which is consistent with the observed trend in the velocity 
difference. In order to explore the link between velocity dispersion and coalescence, we have 
computed 50% probability intervals for the conditional velocity. These are defined to be the 
values of v for which the conditional velocity PDF f{v\r) is fifty percent of its peak value. 
Note that in the absence of velocity dispersion f{v\r) is a delta function centered at {uz\r), 
so the width of the intervals is a measure of dispersion. From the plot of conditional velocity, 
we can note that for large droplets the velocity dispersion is significant. We can also note that 
using DQMOM essentially results in points along the curve {uz\r), i.e., increasing with 
the same choice of moments does not capture the velocity dispersion. 
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Figure 10. Monomodal case with coalescence but no evaporation. Top left: mass density. Top right: 
error on the mass density. Bottom left: velocity difference. Bottom right: error on the velocity differ- 
ence. 



5 Conclusions 



In this work, we have implemented DQMOM to treat the Williams spray equation that de- 
scribes evaporation, acceleration and coalescence of liquid droplets in a laminar gas flow. The 
derivation of the DQMOM equations was shown to be a straightforward task, and resulted in 
a linear system for the source terms. The right-hand side of this linear system is non-zero only 
in the presence of coalescence or non-linear evaporation. The coefficient matrix depends on 
the choice of moments used in DQMOM. 

We have compared this method, as well as the solution obtained with another Eulerian method: 
the multi-fluid model, to the reference solution produced by a classical Lagrangian solver. As 
far as coalescence phenomena are concerned, the efficiency of DQMOM has been shown to 
be better than the multi-fluid model due to its limited numerical diffusion in the size phase 
space, especially for the bimodal distribution function. However, as far as the evaporation 
process is concerned, it is comparable to the multi-fluid model, but still needs a further study 
in order to fully understand how to treat optimally the issue of the evaporative flux due to 
droplet disappearance. Although this issue has been so far neglected in the literature on mo- 
ment methods, our study illustrates that it has an important effect on the moment dynamics. 

The principal conclusion from this study is that DQMOM is numerically robust and straight- 
forward to implement for the Williams spray equation and that it will be a very good candidate 
for more complex two-phase combustion applications once the issue of the evaporative flux 
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is further improved. 
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